1. Introductie: Waarom Statistiek?

1 Introductie

• We leven in het big data era
• Gegevens over locatie, surfgedrag, aankoopgedrag, sociale media
• Wetenschap: expressie meten van duizenden genen, eiwitten,… voor elk subject
• Chemische proces controle: groot aantal sensoren die continu een proces meten zodat het proces kan worden bijgestuurd
• Data driven journalism
• …

Statistiek is de wetenschap om te leren uit empirische gegevens

Statistische geletterdheid is cruciaal om resultaten en publicaties kritisch te kunnen interpreteren.

1.1 Boutade: Met statistiek kan je alles bewijzen

---

1.2 Met foute statistiek kan je alles bewijzen

---

Garcia-Berthou & Alcaraz (Med. Res. Meth. 2004):

• Jaargang 2001 van Nature en British Medical Journal
• 38% en 25% van artikelen bevat een statistische fout
• 11%: fouten tegen interpretatie
• 4%: besluit spreekt evidentie tegen

---

“the expression for many sets of genes was found to be more similar in different tissues within the same species than between species.”

“Overall, our results indicate that there is considerable RNA expression diversity between humans and mice, well beyond what was described previously, likely reflecting the fundamental physiological differences between these two organisms.”

---

https://twitter.com/Y_Gilad/status/593088451462963202

---

https://twitter.com/Y_Gilad/status/593088451462963202

---

2 Oksel microbiome voorbeeld

• https://www.vrt.be/vrtnws/nl/2018/10/22/gezocht-mensen-met-penetrante-lijfgeur-om-probiotische-deodor/
• https://youtu.be/9RIFyqLXdVw

---

• Okselgeur wordt niet veroorzaakt door het zweet zelf. De geur is afkomstig van specificiek micro-organismen van de groep Corynebacterium spp. die het zweet metaboliseren. Staphylococcus spp. zijn een andere groep bacteriën die ook abundant zijn in het microbiome van de oksel die zweet metaboliseren naar metabolieten die niet stinken.

• De CMET-groep aan UGENT doet onderzoek naar het transplanteren van de microbiële gemeenschap, het microbiome, om mensen van geurende oksels af te helpen.

• Voorgestelde therapie
  1. Verwijder het oksel microbiome met antibiotica
  2. Beïnvloed het oksel microbiome door microbiome te transplanteren van een individue die geen geurende oksels heeft (https://youtu.be/9RIFyqLXdVw)

---

2.1 Proefopzet (experimental design)

---

---

• Experiment:

  • 20 personen worden at random uit de populatie getrokken van personen met een okselgeur: steekproef representatief voor populatie!

  • de personen worden at random verdeeld over twee behandelingsgroepen:

    • placebo (enkel antibiotica)
    • transplantie (antibiotica en microbiële transplantatie).
    • Randomisatie is belangrijk om ervoor te zorgen dat de groepen vergelijkbaar zijn.
  • Het microbiome wordt bemonsterd 6 weken na de behandeling.

  • The relative abundantie van Staphylococcus spp. op Corynebacterium spp. + Staphylococcus spp. in het microbiome wordt gemeten via DGGE (Denaturing Gradient Gel Electrophoresis).

---

DGGE

https://doi.org/10.1371/journal.pone.0070538

---

Vertaal onderzoeksvraag naar iets wat we kunnen quantificeren: Is er een verschil in relatieve abundantie van Staphylococcus spp. in het microbiome van de placebo groep en de transplantatie

---

---

2.2 Data Exploratie en beschrijvende statistiek

• Data exploratie is heel belangrijk om inzicht te krijgen in de data en is een essentiële eerste stap om te leren uit data.
• Het wordt vaak ondergewaardeerd of over het hoofd gezien.

---

2.2.1 Importeer de data

• Data in deze cursus wordt verwerkt via het statistisch software pakket R.
• Dit pakket laat toe om te leren uit data.
• In deze cursus gaan we data eerst exploreren om inzicht te verwerven in de gegevens om die vervolgens statistisch te verwerken.

• Vooraleer we hiermee van start kunnen gaan moeten we de data eerst importeren in R.

• Via het volgende commando kunnen we enkele regels van een data bestand inlezen om de structuur van het data bestand te weten te komen.

read_lines("https://raw.githubusercontent.com/statOmics/sbc20/master/data/armpit.csv")

 [1] "trt,rel"                       "placebo,54.99207606973059"    
 [3] "placebo,31.84466019417476"     "placebo,41.09947643979057"    
 [5] "placebo,59.52063914780293"     "placebo,63.573407202216075"   
 [7] "placebo,41.48648648648649"     "placebo,30.44041450777202"    
 [9] "placebo,42.95676429567643"     "placebo,41.7391304347826"     
[11] "placebo,33.896515311510036"    "transplant,57.218124341412015"
[13] "transplant,72.50900360144058"  "transplant,61.89258312020461" 
[15] "transplant,56.690140845070424" "transplant,76"                
[17] "transplant,71.7357910906298"   "transplant,57.757296466973884"
[19] "transplant,65.1219512195122"   "transplant,67.53424657534246" 
[21] "transplant,77.55359394703657" 

• Gegevens in het bestand zijn door comma’s gescheiden.
• Elke rij bevat de gegevens voor 1 proefpersoon
• Verschillende variabelen worden gemeten per persoon en zijn van elkaar gescheiden door een comma. Het bestand is csv formaat: “comma separated values”.
• We kunnen bestanden met dit formaat inlezen R via het commando read_csv.
• We slaan de data op in R in het object met naam ap. Hiervoor gebruiken we de <- operator.
• We geven de data tabel terug door het object aan te roepen door zijn naam te typen.

ap <- read_csv("https://raw.githubusercontent.com/statOmics/sbc20/master/data/armpit.csv")
ap

# A tibble: 20 x 2
   trt          rel
   <chr>      <dbl>
 1 placebo     55.0
 2 placebo     31.8
 3 placebo     41.1
 4 placebo     59.5
 5 placebo     63.6
 6 placebo     41.5
 7 placebo     30.4
 8 placebo     43.0
 9 placebo     41.7
10 placebo     33.9
11 transplant  57.2
12 transplant  72.5
13 transplant  61.9
14 transplant  56.7
15 transplant  76  
16 transplant  71.7
17 transplant  57.8
18 transplant  65.1
19 transplant  67.5
20 transplant  77.6

---

2.2.2 Beschrijvende statistiek

• In artikels en de media worden resultaten uit een steekproef vaak gerapporteerd a.d.h.v. gemiddelde en de standaardafwijking

• We vatten de data eerst samen. We berekenen het gemiddelde en de standaard deviatie (een maat voor de spreiding, zie volgende hoofdstukken). We slaan het resultaat hiervan op in het object apRelSum via apRelSum <-.

1. We pipen (via %>%) het ap dataframe naar de group_by functie om de data te groeperen per treatment trt: group_by(trt).

2. We pipen het resultaat naar de summarize_at function om de “rel” variable samen te vatten en berekenen hierbij het gemiddelde en standaardafwijking. Omdat we de data eerst hebben gegroepeerd zullen we het gemiddelde en de standaard deviatie berekenen per groep.

apRelSum <- ap %>%
  group_by(trt) %>%
  summarize_at("rel",
               list(mean=mean,
                    sd=sd
                    )
                )

We tonen vervolgens het resultaat door het object apRelSum aan te roepen

apRelSum

# A tibble: 2 x 3
  trt         mean    sd
  <chr>      <dbl> <dbl>
1 placebo     44.2 11.5 
2 transplant  66.4  7.88

We kunnen ook een tabel in de webpagina of het pdf bestand integreren via het commando kable van het knitr pakket:

knitr::kable(apRelSum)

trt	mean	sd
placebo	44.15496	11.543251
transplant	66.40127	7.880175

---

2.2.3 Grafieken

We maken in deze cursus gebruik van het pakket ggplot2 om grafieken te maken.
Met de ggplot2 bibliotheek kunnen we gemakkelijk grafieken opbouwen in lagen (layers). Hierdoor leest de code veel makkelijker.

2.2.3.1 barplot

Bar plots worden heel veel gebruikt in artikels om resultaten weer te geven.

1. We pipen de samengevatte data naar de functie ggplot. Dat is de basis van elke ggplot. We selecteren de variabele met de behandeling trt als x variabele en de variabele met naam mean als y-variabele voor de plot. We doen dit steeds via de aestetics aes functie. aes(x=trt,y=mean)

2. We maken een barplot door een laag toe te voegen via de geom_bar function. De statistiek is stat="identity" omdat de hoogte van de bar gelijk is aan de waarde voor y (hier het gemiddelde voor de relatieve abundantie).

3. We voegen foutenvlaggen toe om de onzekerheid op het gemiddelde weer te geven. We doen dit via de geom_errorbar functie en specifiëren het minimum en maximum van de error bar. Het width argument wordt gebruikt om de breedte van de error bar smaller te maken dat deze van de bar.

apRelSum %>%
  ggplot(aes(x=trt,y=mean)) +
  geom_bar(stat="identity") +
  geom_errorbar(aes(ymin=mean-sd,ymax=mean+sd),width=.2)

• Is deze plot informatief??

---

2.2.3.2 boxplots

Barplots zijn geen goede grafieken:

• Ze zijn niet informatief

• Ze visulariseren een samenvatting van de data door twee punten: two point summary: gemiddelde en de standaard deviatie hierop en kunnen beter in een tabel worden opgenomen.

• Ze laten niet toe om anomalieën in de data zoals meetfouten op te sporen.

• Ze verbruiken veel ruimte, b.v. van nul tot de kleinste waarde voor de relatieve abundantie) waar geen data in ligt.

Het is beter om de data zo ruw mogelijk voor te stellen zodat we inzicht krijgen in de verdeling van de gegevens. Hiervoor zullen we ondermeer boxplots gebruiken.

---

2.2.3.3 Boxplots

---

We maken nu een boxplot voor de ap data

1. We pipen het ap dataframe naar ggplot
2. We selecteren de data voor de plot via ggplot(aes(x=trt,y=rel))
3. We voegen laag toe voor de boxplot dmv de functie geom_boxplot()

ap %>%  
  ggplot(aes(x=trt,y=rel)) +
  geom_boxplot()

---

Merk op dat we de plot ook op kunnen slaan als een object.

apBoxplot <- ap %>%  
  ggplot(aes(x=trt,y=rel)) +
  geom_boxplot()

De plot wordt dan niet gemaakt.

Om de plot weer te geven kunnen we het object aanroepen

apBoxplot

---

• Merk op dat we niet zoveel gegevens hebben. Enkel 10 per groep.

• Het is altijd beter om de data zo ruw mogelijk te tonen!

Omdat er niet zoveel gegevens zijn kunnen we de data toevoegen aan de plot zonder dat die te druk wordt.

• Merk op dat we het argument outlier.shape op NA (not available) zetten outlier.shape=NA in the geom_boxplot functie omdat we anders outliers twee keer weer zullen geven. Eerst via de boxplot laag en daarna omdat we een laag met alle ruwe data toevoegen aan de plot.
• We geven de ruwe data weer via de geom_point(position="jitter") functie. We gebruiken hierbij het argument position=‘jitter’ zodat we wat random ruis toevoegen aan de x-cordinaat zodat de gegevens elkaar niet overlappen.

ap %>%  
  ggplot(aes(x=trt,y=rel)) +
  geom_boxplot(outlier.shape=NA) +
  geom_point(position="jitter")

Dit is een informatieve plot!

---

Aangezien we de plot ook hebben opgeslagen konden we dit ook doen door de plot op te roepen en nog een laag toe te voegen.

apBoxplot +
  geom_point(position="jitter")

We hadden het resultaat ook opnieuw op kunnen slaan voor later hergebruik.

apBoxplot <- apBoxplot +
  geom_point(position="jitter")

---

• We zagen duidelijk een effect van de transplantatie op de relatieve abundantie van Staphylococcus.

• Is dat effect nu groot genoeg om te kunnen concluderen dat de behandeling werkt?

---

---

• Inductie: Door middel van statistische besluitvorming (inference) kunnen we uitspraken doen over de populatie op basis van een steekproef.

• De prijs die we hiervoor betalen is onzekerheid!

• We kunnen op basis van een steekproef nooit absoluut zeker zijn van onze conclusies.

---

• Met data kunnen we niet bewijzen dat een behandeling werkt.

• Falsificatie principe van Popper: Data kunnen enkel een hypothese of een theorie ontkrachten.

• Met statistiek kunnen we dus niet aantonen dat de behandeling werkt.

• Statistiek zal ons wel toelaten om het omgekeerde te falsifiëren: als we veronderstellen dat er geen effect van de behandeling, spreekt de data in de steekproef dit tegen?

• Met statistiek kunnen we berekenen hoe waarschijnlijk het is om een gemiddeld verschil in relatieve abundantie te zien tussen placebo en transplantatiegroep dat minstens zo groot is als in de steekproef als de behandeling geen effect zou hebben.

• Die kans wordt een p-waarde genoemd.

• Als p heel klein is, dan is het heel onwaarschijnlijk om een dergelijk effect te observeren in een steekproef door toeval.

• p wordt meestal vergeleken met 5%. Als er geen effect is van de behandeling dan tolereren 5% vals positieve conclusies.

• Om de kans p te berekenen is het nodig om de data te modelleren met statistische modellen.

---

In latere hoofdstukken zullen we zien dat we t-test kunnen gebruiken om hetgeen we observeren in de microbiome dataset te veralgemenen naar de populatie.

ap %>% t.test(rel~trt,data=.)

    Welch Two Sample t-test

data:  rel by trt
t = -5.0334, df = 15.892, p-value = 0.0001249
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -31.62100 -12.87163
sample estimates:
   mean in group placebo mean in group transplant 
                44.15496                 66.40127 

Conclusie: Gemiddeld is de relatieve abundantie van Staphylococcus in het microbiome van personen met een zweetgeur 22.2% hoger na de transplantie dan na de placebo behandeling.

---

2.3 Randomisatie

• Wat wordt zijn de consequenties van het gebruik van een steekproef en van randomisatie?

• Randomisatie is sterk gerelateerd met het concept van de populatie en scope van de studie.

• De scope van de studie moet goed worden omschreven voor de start van het experiment.

• Omdat de statistische analyse valide zou zijn is het noodzakelijk dat de subjecten volledig random worden getrokken uit de populatie naar waar we onze conclusies wensen te veralgemenen.

• Volledig random trekken van de populatie impliceert dat:

  • alle subjecten van de populatie evenveel kans hebben om in de steekproef te worden opgenomen
  • de selectie van een subject onafhankelijk is van de andere subjecten in de steekproef.

• De steekproef is dan representatief voor de populatie, maar is nog steeds random.

• Wat betekent dit?

---

3 Variabiliteit van steekproef tot steekproef

• Om te begrijpen dat een steekproef random is zouden we hetzelfde experiment veel keer moeten kunnen herhalen (repeated sampling).

• Dan zouden we inzicht kunnen krijgen hoe de gegevens veranderen van steekproef tot steekproef.

• Om dit te illustreren zullen we gebruik maken van een hele grote studie.

• Uit die studie zullen we dan herhaaldelijk kleine steekproeven trekken om te begrijpen hoe de gegevens en statistieken veranderen van steekproef tot steekproef. Of om met andere woorden na te gaan wat de variabiliteit is tussen steekproeven.

---

National Health NHanes study

• Sinds 1960 worden elk jaar mensen van alle leeftijden geïnterviewd bij hen thuis.
• Er maakt ook een gezondheidsonderzoek deel uit van de study die in een mobiel onderzoekscentrum wordt afgenomen.
• We zullen deze grote studie gebruiken om at random personen te selecteren van de Amerikaanse populatie.
• Dat zal inzicht geven in hoe de gegevens en resultaten van een analyse zullen variëren van steekproef tot steekproef.
• De data van deze studie is terug te vinden in het R pakket NHANES

---

library(NHANES)
head(NHANES)

# A tibble: 6 x 77
     ID SurveyYr Gender   Age AgeDecade AgeMonths Race1 Race3 Education
  <int> <fct>    <fct>  <int> <fct>         <int> <fct> <fct> <fct>    
1 51624 2009_10  male      34 " 30-39"        409 White <NA>  High Sch…
2 51624 2009_10  male      34 " 30-39"        409 White <NA>  High Sch…
3 51624 2009_10  male      34 " 30-39"        409 White <NA>  High Sch…
4 51625 2009_10  male       4 " 0-9"           49 Other <NA>  <NA>     
5 51630 2009_10  female    49 " 40-49"        596 White <NA>  Some Col…
6 51638 2009_10  male       9 " 0-9"          115 White <NA>  <NA>     
# … with 68 more variables: MaritalStatus <fct>, HHIncome <fct>,
#   HHIncomeMid <int>, Poverty <dbl>, HomeRooms <int>, HomeOwn <fct>,
#   Work <fct>, Weight <dbl>, Length <dbl>, HeadCirc <dbl>, Height <dbl>,
#   BMI <dbl>, BMICatUnder20yrs <fct>, BMI_WHO <fct>, Pulse <int>,
#   BPSysAve <int>, BPDiaAve <int>, BPSys1 <int>, BPDia1 <int>,
#   BPSys2 <int>, BPDia2 <int>, BPSys3 <int>, BPDia3 <int>,
#   Testosterone <dbl>, DirectChol <dbl>, TotChol <dbl>, UrineVol1 <int>,
#   UrineFlow1 <dbl>, UrineVol2 <int>, UrineFlow2 <dbl>, Diabetes <fct>,
#   DiabetesAge <int>, HealthGen <fct>, DaysPhysHlthBad <int>,
#   DaysMentHlthBad <int>, LittleInterest <fct>, Depressed <fct>,
#   nPregnancies <int>, nBabies <int>, Age1stBaby <int>,
#   SleepHrsNight <int>, SleepTrouble <fct>, PhysActive <fct>,
#   PhysActiveDays <int>, TVHrsDay <fct>, CompHrsDay <fct>,
#   TVHrsDayChild <int>, CompHrsDayChild <int>, Alcohol12PlusYr <fct>,
#   AlcoholDay <int>, AlcoholYear <int>, SmokeNow <fct>, Smoke100 <fct>,
#   Smoke100n <fct>, SmokeAge <int>, Marijuana <fct>, AgeFirstMarij <int>,
#   RegularMarij <fct>, AgeRegMarij <int>, HardDrugs <fct>, SexEver <fct>,
#   SexAge <int>, SexNumPartnLife <int>, SexNumPartYear <int>,
#   SameSex <fct>, SexOrientation <fct>, PregnantNow <fct>, gender <dbl>

glimpse(NHANES)

Observations: 10,000
Variables: 77
$ ID               <int> 51624, 51624, 51624, 51625, 51630, 51638, 51646…
$ SurveyYr         <fct> 2009_10, 2009_10, 2009_10, 2009_10, 2009_10, 20…
$ Gender           <fct> male, male, male, male, female, male, male, fem…
$ Age              <int> 34, 34, 34, 4, 49, 9, 8, 45, 45, 45, 66, 58, 54…
$ AgeDecade        <fct>  30-39,  30-39,  30-39,  0-9,  40-49,  0-9,  0-…
$ AgeMonths        <int> 409, 409, 409, 49, 596, 115, 101, 541, 541, 541…
$ Race1            <fct> White, White, White, Other, White, White, White…
$ Race3            <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ Education        <fct> High School, High School, High School, NA, Some…
$ MaritalStatus    <fct> Married, Married, Married, NA, LivePartner, NA,…
$ HHIncome         <fct> 25000-34999, 25000-34999, 25000-34999, 20000-24…
$ HHIncomeMid      <int> 30000, 30000, 30000, 22500, 40000, 87500, 60000…
$ Poverty          <dbl> 1.36, 1.36, 1.36, 1.07, 1.91, 1.84, 2.33, 5.00,…
$ HomeRooms        <int> 6, 6, 6, 9, 5, 6, 7, 6, 6, 6, 5, 10, 6, 10, 10,…
$ HomeOwn          <fct> Own, Own, Own, Own, Rent, Rent, Own, Own, Own, …
$ Work             <fct> NotWorking, NotWorking, NotWorking, NA, NotWork…
$ Weight           <dbl> 87.4, 87.4, 87.4, 17.0, 86.7, 29.8, 35.2, 75.7,…
$ Length           <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ HeadCirc         <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ Height           <dbl> 164.7, 164.7, 164.7, 105.4, 168.4, 133.1, 130.6…
$ BMI              <dbl> 32.22, 32.22, 32.22, 15.30, 30.57, 16.82, 20.64…
$ BMICatUnder20yrs <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ BMI_WHO          <fct> 30.0_plus, 30.0_plus, 30.0_plus, 12.0_18.5, 30.…
$ Pulse            <int> 70, 70, 70, NA, 86, 82, 72, 62, 62, 62, 60, 62,…
$ BPSysAve         <int> 113, 113, 113, NA, 112, 86, 107, 118, 118, 118,…
$ BPDiaAve         <int> 85, 85, 85, NA, 75, 47, 37, 64, 64, 64, 63, 74,…
$ BPSys1           <int> 114, 114, 114, NA, 118, 84, 114, 106, 106, 106,…
$ BPDia1           <int> 88, 88, 88, NA, 82, 50, 46, 62, 62, 62, 64, 76,…
$ BPSys2           <int> 114, 114, 114, NA, 108, 84, 108, 118, 118, 118,…
$ BPDia2           <int> 88, 88, 88, NA, 74, 50, 36, 68, 68, 68, 62, 72,…
$ BPSys3           <int> 112, 112, 112, NA, 116, 88, 106, 118, 118, 118,…
$ BPDia3           <int> 82, 82, 82, NA, 76, 44, 38, 60, 60, 60, 64, 76,…
$ Testosterone     <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ DirectChol       <dbl> 1.29, 1.29, 1.29, NA, 1.16, 1.34, 1.55, 2.12, 2…
$ TotChol          <dbl> 3.49, 3.49, 3.49, NA, 6.70, 4.86, 4.09, 5.82, 5…
$ UrineVol1        <int> 352, 352, 352, NA, 77, 123, 238, 106, 106, 106,…
$ UrineFlow1       <dbl> NA, NA, NA, NA, 0.094, 1.538, 1.322, 1.116, 1.1…
$ UrineVol2        <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ UrineFlow2       <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ Diabetes         <fct> No, No, No, No, No, No, No, No, No, No, No, No,…
$ DiabetesAge      <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ HealthGen        <fct> Good, Good, Good, NA, Good, NA, NA, Vgood, Vgoo…
$ DaysPhysHlthBad  <int> 0, 0, 0, NA, 0, NA, NA, 0, 0, 0, 10, 0, 4, NA, …
$ DaysMentHlthBad  <int> 15, 15, 15, NA, 10, NA, NA, 3, 3, 3, 0, 0, 0, N…
$ LittleInterest   <fct> Most, Most, Most, NA, Several, NA, NA, None, No…
$ Depressed        <fct> Several, Several, Several, NA, Several, NA, NA,…
$ nPregnancies     <int> NA, NA, NA, NA, 2, NA, NA, 1, 1, 1, NA, NA, NA,…
$ nBabies          <int> NA, NA, NA, NA, 2, NA, NA, NA, NA, NA, NA, NA, …
$ Age1stBaby       <int> NA, NA, NA, NA, 27, NA, NA, NA, NA, NA, NA, NA,…
$ SleepHrsNight    <int> 4, 4, 4, NA, 8, NA, NA, 8, 8, 8, 7, 5, 4, NA, 5…
$ SleepTrouble     <fct> Yes, Yes, Yes, NA, Yes, NA, NA, No, No, No, No,…
$ PhysActive       <fct> No, No, No, NA, No, NA, NA, Yes, Yes, Yes, Yes,…
$ PhysActiveDays   <int> NA, NA, NA, NA, NA, NA, NA, 5, 5, 5, 7, 5, 1, N…
$ TVHrsDay         <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ CompHrsDay       <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ TVHrsDayChild    <int> NA, NA, NA, 4, NA, 5, 1, NA, NA, NA, NA, NA, NA…
$ CompHrsDayChild  <int> NA, NA, NA, 1, NA, 0, 6, NA, NA, NA, NA, NA, NA…
$ Alcohol12PlusYr  <fct> Yes, Yes, Yes, NA, Yes, NA, NA, Yes, Yes, Yes, …
$ AlcoholDay       <int> NA, NA, NA, NA, 2, NA, NA, 3, 3, 3, 1, 2, 6, NA…
$ AlcoholYear      <int> 0, 0, 0, NA, 20, NA, NA, 52, 52, 52, 100, 104, …
$ SmokeNow         <fct> No, No, No, NA, Yes, NA, NA, NA, NA, NA, No, NA…
$ Smoke100         <fct> Yes, Yes, Yes, NA, Yes, NA, NA, No, No, No, Yes…
$ Smoke100n        <fct> Smoker, Smoker, Smoker, NA, Smoker, NA, NA, Non…
$ SmokeAge         <int> 18, 18, 18, NA, 38, NA, NA, NA, NA, NA, 13, NA,…
$ Marijuana        <fct> Yes, Yes, Yes, NA, Yes, NA, NA, Yes, Yes, Yes, …
$ AgeFirstMarij    <int> 17, 17, 17, NA, 18, NA, NA, 13, 13, 13, NA, 19,…
$ RegularMarij     <fct> No, No, No, NA, No, NA, NA, No, No, No, NA, Yes…
$ AgeRegMarij      <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 20,…
$ HardDrugs        <fct> Yes, Yes, Yes, NA, Yes, NA, NA, No, No, No, No,…
$ SexEver          <fct> Yes, Yes, Yes, NA, Yes, NA, NA, Yes, Yes, Yes, …
$ SexAge           <int> 16, 16, 16, NA, 12, NA, NA, 13, 13, 13, 17, 22,…
$ SexNumPartnLife  <int> 8, 8, 8, NA, 10, NA, NA, 20, 20, 20, 15, 7, 100…
$ SexNumPartYear   <int> 1, 1, 1, NA, 1, NA, NA, 0, 0, 0, NA, 1, 1, NA, …
$ SameSex          <fct> No, No, No, NA, Yes, NA, NA, Yes, Yes, Yes, No,…
$ SexOrientation   <fct> Heterosexual, Heterosexual, Heterosexual, NA, H…
$ PregnantNow      <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ gender           <dbl> 1, 1, 1, 1, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1,…

---

3.1 Data exploratie

Onderzoeksvraag: hoe verschilt de lengte van volwassen mannen en vrouwen.

1. We pipen de dataset naar de function filter om de data te filteren volgens leeftijd.
   
2. We plotten de lengte metingen.
   • We selecteren de data data met het commando ggplot(aes(x=lengte))
   • We voegen een histogram toe met het commando geom_histogram()
   • We maken twee vertikale panels met het commando facet_grid(Gender~.)
   • We veranderen het label van de x-as met de xlab functie.

NHANES%>%  
  filter(Age > 18 & !is.na(Height)) %>%
  ggplot(aes(x = Height))+
  geom_histogram() +
  facet_grid(Gender ~ .) +
  xlab("Lengte (cm)")

---

We zien dat de data nu min of meer symmetrisch verdeeld zijn en een klokvorm hebben.
Dat zal ons toe laten om de data verder samen te vatten door gebruik te maken van twee statistieken: het gemiddelde en de standaard deviatie wat een maat is voor de spreiding van de gegevens rond het gemiddelde.

---

We maken nu een subset van de data die we zullen gebruiken om aan te tonen hoe de variabiliteit in kleine steekproeven kan variëren van steekproef tot steekproef.

1. We filteren op leeftijd en verwijderen ontbrekenden gegevens (NA, Not Available).
2. We selecteren enkel het geslacht en Lengte zodat de dataset geen onnodige variabelen bevat.

nhanesSub <- NHANES %>%
  filter(Age > 25 & !is.na(Height)) %>%
  select(c("Gender","Height"))

---

We berekenen het gemiddelde en de standaard deviatie voor de lengte voor mannen en vrouwen in de grote dataset. We groeperen de data hiervoor op basis van het geslacht (variable Gender).

HeightSum <- nhanesSub %>%
  group_by(Gender) %>%
  summarize_at("Height",
               list(mean = mean,
               sd = sd)
              )

knitr::kable(
  HeightSum %>%
  mutate_if(is.numeric, round, digits=1)
  )

Gender	mean	sd
female	161.9	7.2
male	175.8	7.5

---

3.2 Experiment

• Stel dat we geen toegang hebben tot de metingen van de NHANES studie.

• We zouden dan een experiment op moeten zetten om metingen bij mannen en vrouwen te doen.

• Veronderstel dat we budget hebben om metingen bij 5 mannen en 5 vrouwen te doen.

• We zouden dan 5 mannen en 5 vrouwen boven de 25 jaar at random selecteren uit de Amerikaanse populatie.

• We kunnen dit experiment simuleren door 5 vrouwen en 5 mannen at random te selecteren uit de NHANES studie.

---

set.seed(1000001)
nSamp <- 5
fem <- nhanesSub %>%
  filter(Gender=="female") %>%
  sample_n(size=5)

mal <- nhanesSub %>%
  filter(Gender=="male") %>%
  sample_n(size=5)

samp1 <- rbind(fem,mal)

samp1

# A tibble: 10 x 2
   Gender Height
   <fct>   <dbl>
 1 female   159.
 2 female   168.
 3 female   157.
 4 female   155.
 5 female   163.
 6 male     178.
 7 male     173.
 8 male     187.
 9 male     173.
10 male     185.

---

Data Exploratie

samp1 %>%
  ggplot(aes(x=Height)) +
  geom_histogram() +
  facet_grid(Gender~.) +
  xlab("Lengte (cm)")

HeightSumExp1 <- samp1 %>%
  group_by(Gender) %>%
  summarize_at("Height",
               list(mean = mean,
                      sd = sd)
                  )
HeightSumExp1

# A tibble: 2 x 3
  Gender  mean    sd
  <fct>  <dbl> <dbl>
1 female  161.  5.17
2 male    179.  6.49

Histogram is niet zinvol als we maar zo weinig datapunten hebben.

---

Boxplot is beter:

---

We voeren hier ook een t-test uit.

t.test(Height~Gender,data=samp1)

    Welch Two Sample t-test

data:  Height by Gender
t = -4.9713, df = 7.6177, p-value = 0.001261
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -27.068985  -9.811015
sample estimates:
mean in group female   mean in group male 
              160.62               179.06 

In het experiment zijn vrouwen zijn gemiddeld 18.44 cm kleiner dan mannen. En als we een statistische test uitvoeren (zie hoofdstuk 5: Statistische besluitvorming) kunnen we besluiten dat dit verschil statistisch significant is.

---

3.3 Herhaal het experiment

Als we het experiment herhalen selecteren we andere mensen en verkrijgen we andere resultaten.

fem <- nhanesSub %>%
  filter(Gender=="female") %>%
  sample_n(size=5)

mal <- nhanesSub %>%
  filter(Gender=="male") %>%
  sample_n(size=5)

samp2 <- rbind(fem,mal)

HeightSumExp2 <- samp2 %>%
  group_by(Gender) %>%
  summarize_at("Height",
               list(mean=mean,
                    sd=sd)
              )
HeightSumExp2

# A tibble: 2 x 3
  Gender  mean    sd
  <fct>  <dbl> <dbl>
1 female  169.  9.38
2 male    175.  6.12

samp2 %>%
  ggplot(aes(x = Gender,y = Height)) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = "jitter") +
  geom_point(
    aes(x = 1, y = HeightSumExp1$mean[1]),
    size = 3,
    pch = 17,
    color="darkred") +
  geom_point(
    aes(x = 2, y = HeightSumExp1$mean[1]),
    size = 3,
    pch = 17,
    color = "darkred") +
  ylab("Height (cm)")

t.test(Height ~ Gender, data=samp2)

    Welch Two Sample t-test

data:  Height by Gender
t = -1.1379, df = 6.8813, p-value = 0.2932
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -17.586435   6.186435
sample estimates:
mean in group female   mean in group male 
              169.06               174.76 

In de nieuwe steekproef zijn vrouwen zijn gemiddeld 5.7 cm kleiner dan mannen. En dit verschil is statistisch niet significant

---

3.4 Herhaal het experiment opnieuw

seed <- 1336602
fem <- nhanesSub %>%
  filter(Gender=="female") %>%
  sample_n(size=5)

mal <- nhanesSub %>%
  filter(Gender=="male") %>%
  sample_n(size=5)

samp3 <- rbind(fem,mal)

HeightSumExp3 <- samp3 %>%
  group_by(Gender) %>%
  summarize_at("Height",
               list(mean=mean,
                    sd=sd)
              )
HeightSumExp3

# A tibble: 2 x 3
  Gender  mean    sd
  <fct>  <dbl> <dbl>
1 female  162.  12.5
2 male    170.  14.5

samp3 %>%
  ggplot(aes(x = Gender,y = Height)) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = "jitter") +
  geom_point(
    aes(x = 1, y = HeightSumExp1$mean[1]),
    size = 3,
    pch = 17,
    color="darkred") +
  geom_point(
    aes(x = 2, y = HeightSumExp1$mean[1]),
    size = 3,
    pch = 17,
    color = "darkred") +
  ylab("Height (cm)")

t.test(Height ~ Gender, data=samp3)

    Welch Two Sample t-test

data:  Height by Gender
t = -0.99138, df = 7.8343, p-value = 0.3511
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -28.34441  11.34441
sample estimates:
mean in group female   mean in group male 
              161.96               170.46 

In de nieuwe steekproef zijn vrouwen gemiddeld 8.5 cm kleiner dan mannen. En dit verschil is statistisch niet significant

---

3.5 Samenvatting

• We trokken at random andere proefpersonen in elke steekproef

• Hierdoor verschillen lengtemetingen van steekproef tot steekproef.

• Dus ook de geschatte gemiddeldes en standaard deviaties.

• Bijgevolg zijn onze conclusies ook onzeker en kunnen deze wijzigen van steekproef tot steekproef.

• Steekproeven waarbij het effect tegengesteld is aan dat in de populatie en waarbij we besluiten dat het verschil significant is, zijn zeldzaam.

\(\rightarrow\) Met statistiek controleren we de kans op het trekken foute conclusies.

• Dit wordt in de onderstaande code geïllustreerd. We trekken 10000 herhaalde steekproeven van 5 vrouwen en 5 mannen uit de NHANES studie.

# Aantal simulaties en steekproefgrootte per groep
nSim <- 10000
nSamp <- 5

# We filteren de data vooraf zodat we dit niet telkens opnieuw hoeven te doen
fem <- nhanesSub %>%
  filter(Gender == "female")

mal <- nhanesSub %>%
  filter(Gender == "male")

# Simulatie studie
# Om snelle functies te kunnen gebruiken nemen we eerst nSim steekproeven en berekenen we daarna alles.

femSamps <- malSamps <-matrix(NA, nrow=nSamp, ncol=nSim)
for (i in 1:nSim)
{
  femSamps[,i] <- sample(fem$Height, nSamp)
  malSamps[,i] <- sample(mal$Height, nSamp)
}

res <- data.frame(
  verschil=colMeans(femSamps) - colMeans(malSamps),
  Rfast::ttests(femSamps, malSamps)
  )

sum(res$pvalue < 0.05 & res$verschil < 0)

[1] 7148

sum(res$pvalue >= 0.05)

[1] 2852

sum(res$pvalue < 0.05 & res$verschil>0)

[1] 0

res %>%
  ggplot(aes(x=verschil,y=-log10(pvalue),color=pvalue < 0.05)) +
  geom_point() +
  xlab("Gemiddeld Verschil (cm)") +
  ylab("Statistische Significantie (-log10 p)")

res %>%
  ggplot(aes(y = verschil)) +
  geom_boxplot() +
  ylab("Gemiddeld Verschil (cm)")

  xlab("")

$x
[1] ""

attr(,"class")
[1] "labels"

Op basis van 10 000 steekproeven van 5 mannen en 5 vrouwen zagen we dat in 7148 steekproeven vrouwen gemiddeld significant kleiner zijn dan mannen. In 2852 steekproeven besluiten we dat vrouwen en mannen gemiddeld niet significant verschillen in lengte. En in 0 besluiten we dat vrouwen gemiddeld significant groter zijn dan mannen.

• De steekproef die we toonden waaruit we zouden besluiten dat vrouwen significant groter zijn dan mannen is heel onwaarschijnlijk. Er moesten 1336602 steekproeven worden getrokken om deze extreme steekproef te vinden.

Het feit dat we in veel steekproeven resultaten vinden die statistisch niet significant zijn komt omdat de statistische toets een te lage kracht heeft om het verschil te detecteren wanneer er maar 5 observaties zijn per groep.

---

3.6 Grotere steekproef?

Wat gebeurt er als we de steekproef verhogen naar 50 observaties per groep?

# Aantal simulaties en steekproefgrootte per groep
nSim <- 10000
nSamp <- 50

# We filteren de data vooraf zodat we dit niet telkens opnieuw hoeven te doen
fem <- nhanesSub %>%
  filter(Gender == "female")

mal <- nhanesSub %>%
  filter(Gender == "male")

# Simulatie studie
# Om snelle functies te kunnen gebruiken nemen we eerst nSim steekproeven en berekenen we daarna alles.

femSamps <- malSamps <- matrix(NA, nrow = nSamp, ncol = nSim)
for (i in 1:nSim)
{
  femSamps[,i] <- sample(fem$Height, nSamp)
  malSamps[,i] <- sample(mal$Height, nSamp)
}

res <- data.frame(
  verschil = colMeans(femSamps) - colMeans(malSamps),
  Rfast::ttests(femSamps, malSamps)
  )

sum(res$pvalue < 0.05 & res$verschil < 0)

[1] 10000

sum(res$pvalue >= 0.05)

[1] 0

sum(res$pvalue < 0.05 & res$verschil > 0)

[1] 0

res %>%
  ggplot(aes(x=verschil,y=-log10(pvalue),color=pvalue<0.05)) +
  geom_point() +
  xlab("Gemiddeld Verschil (cm)") +
  ylab("Statistische Significantie (-log10 p)")

res %>%
  ggplot(aes(y=verschil)) +
  geom_boxplot() +
  ylab("Gemiddeld Verschil (cm)")

  xlab("")

$x
[1] ""

attr(,"class")
[1] "labels"

---

4 Salk Study

• In 1916, brak de eerste grote polio epidemie uit in de USA.
• Begin de jaren 50 ontwikkelde John Salk een vaccin met belovende resultaten in het lab.
• In 1954, heeft de National Foundation for Infantile Paralysis (NFIP) een grote studie opgezet om de effectiviteit van het Salk vaccin na te gaan.
• Veronderstel dat de NFIP in 1954 een groot aantal kinderen zou hebben gevaccineerd, wat zouden ze dan kunnen besluiten als de polio incidentie in 1954 lager was dan in 1953?

---

4.1 NFIP Study

4.1.1 Design

• Grote simultane studie met gevaccineerde kinderen (cases) en ongevaccineerde kinderen (controles).
• In scholen van districten met hoge polio incidentie.
• Cases: kinderen van de tweede graad van het lager onderwijs waarvan de ouders toestemden met vaccinatie.
• Controles: kinderen van de eerste en derde graad.

4.1.2 Data

nfip <- tibble(
  group=c("cases","controls","noConcent"),
  grade=c("g2","g1g3","g2"),
  vaccin=c("yes","no","no"),
  total=c(221998,725173,123605),
  polio=c(54,391,56)
  ) %>%
  mutate(noPolio = total - polio)
knitr::kable(nfip)

group	grade	vaccin	total	polio	noPolio
cases	g2	yes	221998	54	221944
controls	g1g3	no	725173	391	724782
noConcent	g2	no	123605	56	123549

Vergelijk de polio incidentie?

---

nfip <- nfip %>%
  mutate(incidencePM = round(nfip$polio/nfip$total*1e6,0))
knitr::kable(nfip)

group	grade	vaccin	total	polio	noPolio	incidencePM
cases	g2	yes	221998	54	221944	243
controls	g1g3	no	725173	391	724782	539
noConcent	g2	no	123605	56	123549	453

Wat kunnen we concluderen?

---

4.2 Confounding

• We observeren een lagere polio (P) incidentie voor kinderen bij wie de ouders geen toestemming gaven dan in de controle groep.

• Toestemming voor vaccinatie (V) is geassocieerd met de socio-economische status (S).

• Kinderen van lagere socio-economische status zijn meer resistent tegen de ziekte.

• De groepen van cases en controles zijn niet vergelijkbaar
  • verschil in leeftijd
  • verschil in socio-economische status en
  • verschil in vatbaarheid voor de ziekte.

---

4.3 Salk Study

4.3.1 Design

Een nieuwe studie werd uitgevoerd: dubbel blinde gerandomiseerde studie.

• Kinderen worden at random toegewezen aan controle of case arm van het experiment nadat de ouders toestemden met vaccinatie.
• Controle: vaccinatie met placebo
• Treatment: vaccinatie met vaccin
• Double blinding:
  • ouders en kinderen weten niet of ze werden gevaccineerd of niet
  • medische staf en onderzoekers weten niet of het kind het vaccin of de placebo kreeg

---

4.3.2 Data

salk <- data.frame(
  group=c("cases","control","noConcent"),
  treatment=c("vaccine","placebo","none"),
  total=c(200745,201229, 338778),polio=c(57,142,157)
  ) %>%
  mutate(
    noPolio = total-polio,
    incidencePM = round(polio/total*1e6,0)
    )
knitr::kable(salk)

group	treatment	total	polio	noPolio	incidencePM
cases	vaccine	200745	57	200688	284
control	placebo	201229	142	201087	706
noConcent	none	338778	157	338621	463

• We observeren een veel groter effect nu dat cases en controles vergelijkbaar zijn, incidentie van respectievelijk 284 and 706 per miljoen.

• De polio incidentie voor kinderen die geen toestemming geven blijft vergelijkbaar 453 and 463 per miljoen respectievelijk in the NFIP and Salk study.

---

5 Rol van Statistiek

• We hebben gezien dat
  • het belangrijk is om de scope van de studie goed te specifiëren voor de start van het experiment
  • randomisatie nodig is om een representatieve steekproef te nemen
  • steekproef grootte is heel belangrijk
  • we moeten ons bewust zijn van Confounding
  • een goede controle is belangrijk

\(\rightarrow\) Goede proefopzet is cruciaal!

---

• We hebben ook geobserveerd dat er variabiliteit is in de populatie
• We kunnen maar een beperkte steekproef nemen uit de populatie

\(\rightarrow\) onzekerheid in de resultaten \(\rightarrow\) onzekerheid in de conclusies

---

• Statistiek is de wetenschap voor het
  1. verzamelen (experimenteel design),
  2. exploren (data exploration) en
  3. leren van data zodat we hetgeen we observeren in de steekproef zouden kunnen veralgemenen naar de populatie terwijl we de onzekerheid quantificeren, controleren en rapporteren (statistisch modelleren en statistische besluitvorming).
• Statistiek speelt daarom een heel belangrijke rol in zowat alle wetenschappen

Home